30道练习题带你玩转统计学的R语言版

前天在天津医科大学做生物信息学分享的时候，提到了小伊老师的统计学课程，不过昨天航班晚点，来不及整理分享回答大家的疑惑，所以今天才出这30个题目。

统计学是一门很深的学问，这里仅仅是出题帮助大家熟练使用R语言来学习统计学知识，具体知识点需要更深入阅读书籍或者教程：

• 推荐一下 统计学基础：https://mp.weixin.qq.com/s/OtB2h6f00U2SRZLzveJKfQ

• 统计学精华-statQuest教学视频：https://mp.weixin.qq.com/s/X0PE9S0BgSuCcAV9zeY1jQ

基础概念

需要掌握R内置数据集及R包数据集

• 内置数据集：https://mp.weixin.qq.com/s/dZPbCXccTzuj0KkOL7R31g

• airway 数据集 https://mp.weixin.qq.com/s/FTdOyxvlS4L-zL3inHTtZg

理解 定性变量（qualitative variable） 和 定量变量(quantitative variable)

定量数据的集中趋势指标主要是：众数、分位数和平均数

定量数据的离散趋势指标主要是：极差，方差和标准差，标准分数，相对离散系数（变异系数），偏态系数与峰态系数

Q1: 载入R中自带的数据集 iris，指出其每列是定性还是定量数据

Q2: 对数据集 iris的所有定量数据列计算集中趋势指标：`众数、分位数和平均数`

Q3:对数据集 iris的所有定性数据列计算水平及频次

Q4:对数据集 iris的所有定量数据列计算离散趋势指标：方差和标准差等

Q5:计算数据集 iris的前两列变量的相关性，提示cor函数可以选择3种methods

Q6:对数据集 iris的所有定量数据列内部zcore标准化，并计算标准化后每列的平均值和标准差

Q7:计算列内部zcore标准化后 iris的前两列变量的相关性

Q8: 根据数据集 iris的第五列拆分数据集后重复上面的Q2到Q7问题

Q9:载入R中自带的数据集 mtcars，重复上面的Q1到Q7个问题

Q10: 载入r包airway并且通过assay函数拿到其表达矩阵后计算每列之间的相关性

关于 airway 代码如下，需要理解：

options(stringsAsFactors = F)
library(airway)
data(airway)
# 这里需要自行学习bioconductor里面的RangedSummarizedExperiment对象
airway 
RNAseq_expr=assay(airway)
colnames(RNAseq_expr) 
RNAseq_expr[1:4,1:4] 
# RNAseq_expr 是一个数值型矩阵，属于连续性变量，可以探索众数、分位数和平均数 ，极差，方差和标准差等统计学指标
RNAseq_gl=colData(airway)[,3]
table(RNAseq_gl)

是 8个样本的RNA-seq数据的counts矩阵，这8个样本分成2组，每组是4个样本， 分别是 trt 和 untrt 组。

通过上面的代码，我们得到了对airway数据集的RNA-seq数据的counts矩阵，命名为 RNAseq_expr 下面会用得到。

表达矩阵相关

首先了解各种统计分布：https://mp.weixin.qq.com/s/uly4jlQomk9LZlHyknkNdg 在R语言的实现方式。

Q1: 把RNAseq_expr第一列全部加1后取log2后计算平均值和标准差

tmp=log2(RNAseq_expr[,1]+1)
mean(tmp)
sd(tmp)

Q2: 根据上一步得到平均值和标准差生成同样个数的随机的正态分布数值

a=rnorm(length(tmp),mean = mean(tmp),sd = sd(tmp))
a=sort(a)
plot(a)
points(sort(tmp))

Q3: 删除RNAseq_expr第一列低于5的数据后，重复Q1和Q2

tmp=RNAseq_expr[,1]
tmp=tmp[tmp>5]
tmp=log2(tmp)
a=rnorm(length(tmp),mean = mean(tmp),sd = sd(tmp))
a=sort(a)
plot(a)
points(sort(tmp))

Q4: 基于Q3对RNAseq_expr的第一列和第二列进行T检验

x=RNAseq_expr[,1]
x=x[x>5]
x=log2(x)

y=RNAseq_expr[,2]
y=y[y>5]
y=log2(y)

t.test(x,y)
library(ggpubr)
df=data.frame(value=c(x,y),
              group=c(rep('x',length(x)),rep('y',length(y))))
ggboxplot(df, y = "value", x = "group")

Q5:     取RNAseq_expr行之和最大的那一行根据分组矩阵进行T检验

pos=which.max(rowSums(RNAseq_expr))
t.test(RNAseq_expr[pos,]~RNAseq_gl)
pos

Q6: 取RNAseq_expr的MAD最大的那一行根据分组矩阵进行T检验

pos=which.max(apply(RNAseq_expr,1,mad))
t.test(RNAseq_expr[pos,]~RNAseq_gl)
pos

Q7: 对RNAseq_expr全部加1后取log2后重复Q5和Q6

RNAseq_expr=log2(RNAseq_expr+1)
pos=which.max(rowSums(RNAseq_expr))
pos
t.test(RNAseq_expr[pos,]~RNAseq_gl)
pos=which.max(apply(RNAseq_expr,1,mad))
pos
t.test(RNAseq_expr[pos,]~RNAseq_gl)

看看是不是基因变化了，统计结果也变化了

Q8: 取RNAseq_expr矩阵的MAD最高的100行，对列和行分别进行层次聚类

cg=names(tail(sort(apply(RNAseq_expr,1,mad)),100))
dat=RNAseq_expr[cg,]
plot(hclust(dist(t(dat))))
colnames(dat)
RNAseq_gl
plot(hclust(dist( dat )))

检查一下聚类结果跟样本的处理信息是否能对应

Q9: 取RNAseq_expr矩阵的SD最高的100行，对列和行分别进行层次聚类

Q10: 对Q8矩阵按照行和列分别归一化并且热图可视化

cg=names(tail(sort(apply(RNAseq_expr,1,mad)),100))
dat=RNAseq_expr[cg,]
pheatmap::pheatmap(scale(dat))
pheatmap::pheatmap(t(scale(t(dat))))

统计检验相关

这里需要对前面的RNAseq_expr矩阵进行一定程度的过滤，主要是过滤那些每一列都为0的行。

options(stringsAsFactors = F)
rm(list=ls())
library(airway) 
RNAseq_expr=assay(airway)
e1=RNAseq_expr[apply(RNAseq_expr,1,function(x) sum(x>0)>1),] 
colnames(RNAseq_expr) 
RNAseq_gl=colData(airway)[,3]
table(RNAseq_gl) 

上面的 e1 矩阵下面的习题就用得到。

Q1: 对e1每一行独立根据分组矩阵进行T检验，检查为什么有些行T检验失败

apply(e1, 1, function(x){
  t.test(x~RNAseq_gl)$p.value
}) 

Q2: 找出T检验失败的行并且从e1矩阵剔除掉

e1_a=e1[,RNAseq_gl=='trt']
e1_b=e1[,RNAseq_gl=='untrt']
a_filter=apply(e1_a, 1,function(x) sd(x)>0)
b_filter=apply(e1_b, 1,function(x) sd(x)>0)
table(a_filter,b_filter)
e1=e1[a_filter | b_filter,]

Q3: 对过滤后的e1矩阵进行每一行独立根据分组矩阵进行T检验

Q4: 对e1矩阵进行加1后log2的归一化命名为e2再对每一行独立根据分组矩阵进行T检验

Q5: 对e1,e2的T检验P值做相关性分析

p1=apply(e1, 1, function(x){
  t.test(x~RNAseq_gl)$p.value
}) 
e2=log(e1+1)
p2=apply(e2, 1, function(x){
  t.test(x~RNAseq_gl)$p.value
}) 
plot(p1,p2)
cor(p1,p2)

写到这里，我突然间感觉自己代码很壮观，也不知道为什么，就觉得自己写代码很帅气，后面的5个题目就不写啦，这些已经够大家用的。

我这里并没有提到基因和样本这样的词语，就是希望其他领域学R的朋友也可以看看，如果生物信息学领域，这样的简单T检验是有很多不合理的地方，比如文库大小，比如基因表达量分布等等。

本习题是生信工程师全套教学视频（74小时）的一个补充

不少表现积极的朋友提交了自己的作业笔记和学习心得，非常棒！

和10万人一起学生信

。。。此列表待你更新。。。。